arXiv:1502.02280vl [cs.NA] 8 Feb 2015 


Noname manuscript No. 

(will be inserted by the editor) 


A comparison of the Extrapolated Successive 
Overrelaxation and the Preconditioned Simultaneous 
Displacement methods for augmented linear systems 


M. A. Louka • N. M. Missirlis 


Received: date / Accepted: date 


Abstract In this paper we study the impact of two types of precondition¬ 
ing on the numerical solution of large sparse augmented linear systems. The 
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ative method to converge. In addition, we develop a geometric approach, for 
determining the optimum values of their parameters and corresponding spec¬ 
tral radii. It is shown that both iterative methods studied (GMESOR and 
GMPSD) attain the same rate of convergence. Numerical results confirm our 
theoretical expectations. 
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1 Introduction 

Let A € be a symmetric positive definite matrix and B € be a 

matrix of full column rank, where m > n. Then, the augmented linear system 
is of the form 


Au = b 


( 1 ) 


where 


A = 






( 2 ) 


with denoting the transpose of the matrix B. 

When A and B are large and sparse matrices, iterative methods for solving 
0-0 are effective and more attractive than direct methods, because of stor¬ 
age requirements and preservation of sparsity. There are several approaches 
to the iterative solution of 0-0. First, we mention multigrid methods |13) . 
|35) . which are often the most efficient and effective methods for solving large, 
sparse, linear systems m, EH]- For example, one can apply multigrid tech¬ 
niques to the whole system 0-0 to solve problems in areas of computational 
fluid dynamics [52] , [30] , [10] , [S] , [S3] , [31] , [3S] constrained optimization sa, 
[44] , [45] , [46] , mixed finite elements [I] , [23] and elsewhere. For parallel multi¬ 
grid see e.g [20], [32], ES- 

On the other hand the difficulty in applying iterative methods such as the 
Successive Overrelaxation (SOR) method [5S] to the system 0-0 is the sin¬ 
gularity of the block diagonal part of the coefficient matrix. Various methods 
have been developed to overcome this problem such as the Uzawa and the 
Preconditioned Uzawa methods [2], [H], [IS]- In 2001, Golub et al. [21] gener¬ 
alized the Uzawa and the Preconditioned Uzawa methods by introducing an 
additional acceleration parameter and produced the SOR-like method. When 
a good preconditioning matrix is easily computed one can consider the MIN- 
RES and GMRES methods [2Q], [M] for solving 0-0. In case the matrix A 
in 0 is symmetric and positive definite, the Preconditioned Conjugate Gra¬ 
dient (PCG) method [31] can be applied. This was done with an SOR-like 
preconditioner in the work by Li, Evans and Zhang in [33] . In 2005, Bai et al. 
[5] studied the Generalized SOR (GSOR) method by introducing an additional 
parameter to the SOR-like method and proved that it possesses the same rate 
of convergence but lower complexity than the PCG method. Furthermore, the 
Generalized Modified Extrapolated SOR (GMESOR) method was also pro¬ 
posed for further study. The latter is a generalization of the GSOR method as 
it uses one additional parameter. The way of introducing parameters resembles 
the one followed for the formulation of the Modified SOR method [55], [29] . 
[36] . [37] . [38] in case of two-cyclic linear systems. 

The present paper was motivated by the work in [5]. We develop the con¬ 
vergence analysis of the Generalized Modified Extrapolated SOR (GMESOR) 
method and the Generalized Modified Preconditioned Simultaneous Displace¬ 
ment (GMPSD) method. These methods introduce more parameters with the 
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hope to further increase their rate of convergence. The goal of our work was 
to study the impact of two different preconditioning matrices to the conver¬ 
gence rate of the associated iterative method for solving the augmented linear 
system ©-©• First, we use the preconditioning matrix which is formed by 
the lower triangular part of A and formulate the GMESOR method which 
is an extrapolated form of the GSOR method. Secondly, we consider as pre¬ 
conditioning matrix the product of the lower with the upper triangular part 
of A and construct the GMPSD method. The reason for studying the latter 
form of preconditioning matrix is to obtain a better approximation to the 
matrix A than the former and as such it is hoped to produce an iterative 
method with a faster rate of convergence. The construction of both meth¬ 
ods resembles the one followed for the MESOR and MPSD methods studied 
in [87) and [38], respectively, for two-cyclic matrices. Our starting point, for 
studying these iterative methods, is the derivation of functional relationships 
which relate the eigenvalues of their iteration matrices with those of the ma¬ 
trix J = A~^B. Assuming that the matrix Q is symmetric positive or 

negative definite, the eigenvalues of the matrix J are real and either positive or 
negative, respectively. Under these assumptions we find sufficient conditions 
for the convergence of the GMESOR and GMPSD methods and determine 
the optimum values of their parameters. The study of GMESOR and GMPSD 
becomes interesting as these methods can be used either as preconditioners 
to Krylov subspace methods H, HB, ED, m or as smoothers to multilevel 
methods El’ El’ El- Traditionally, multigrid methods utilize stationary iter¬ 
ative methods (such as Jacobi, Gauss-Seidel ) to smooth out high-frequency 
errors and accelerate the convergence. In [32] a semi-iterative method, namely 
the Chebyshev-Jacobi method, was used as smoother. Similarly, the GMESOR 
method or the GMPSD method in combination with semi-iterative techniques 
can be used as smoothers to accelerate the rate of convergence of multigrid 
methods. Recent work for the application of algebraic multigrids for saddle 
point systems is presented in |35j and the references therein. 

The contributions of our paper can be summarized as follows. 

(i) We present a simple and unified approach for developing the convergence 
analysis of the GMESOR and GMPSD methods. In particular, we develop a 
geometrical approach for the determination of the optimum values of the pa¬ 
rameters in GMESOR and GMPSD methods which is similar to Varga EO] p. 
Ill, for finding the optimum value of the parameter oj in SOR. The difference, 
in our case, is that now the functional relationship contains more than two 
parameters and consequently we had to extend the proof of EOj. There is a 
different algebraic approach in [55] pp. 279 for the determination of the two 
optimum values for uj and w in the Modified SOR (MSOR) method which, 
with some additional modifications, will solve the problem as far as the GSOR 
method is concerned. Nevertheless, it is doubtful whether this approach works 
also for the determination of the optimum value for more than two parameters 
as is the case for the GMESOR and GMPSD methods. This is also the case if 
one adopts the approach of |S] . 

(ii) Prom our theoretical and experimental analysis it is shown that both afore- 
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mentioned forms of preconditioning matrices have the same impact on the 
convergence rate of the induced iterative method for the numerical solution 
of the augmented linear system (CD-©. More specifically, the GMESOR and 
GMPSD methods attain the same convergence rate since their spectral radii 
are identical for the optimum values of their parameters. In particular, we 
show that GMESOR degenerates to GSOR, whereas a simplified version of 
the GMPSD method is identical to a backward form of the GSOR method. 
Furthermore, we compare the effectiveness of our methods in relation to the 
PHSS [5], [7], [S], [in], [H], [Sm and Krylov subspace methods [H], [32], [IS]- 
The paper is organized as follows. In section 2 we study the convergence of the 
GMESOR method. In particular, we find sufficient conditions for GMESOR to 
converge under the assumption that the eigenvalues of the J matrix are real. 
We also determine optimum values for its parameters. A similar convergence 
analysis for the GMPSD method is developed in section 3. In section 4, we 
present our numerical results and finally in section 5 we state our remarks and 
conclusions. 


2 The Generalized Modified Extrapolated SOR (GMESOR) 
method 

Let the coefficient matrix M of (HJ be defined by the splitting 


A=V-C-U 


(3) 


where 



/ 0 0 \ 
\B'^ aQj 



-B \ 

a-a)Qj ’ 


(4) 


with Q G be a prescribed nonsingular and symmetric matrix and a £ K.. 

Furthermore, we denote by T, the diagonal matrix T = diag{TiIm,T 2 ln) with 
Ti, r 2 G R — {0}, Im G and G be identity matrices. 

For the numerical solution of ©, we consider the following iterative scheme 


(yUi)) =^(^ 1 -^ 2 ) +v{ri,T2) (5) 

where 

n{Ti,T2)=I-R-^TA, r,{Ti,T2)=R-^Tb, ( 6 ) 

i? is a nonsingular matrix to be defined and I = diag{Im, In)- 

In the sequel we consider two different types of preconditioning matrices R 

and study the corresponding iterative methods derived by and 
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2.1 The functional relationship 

As a first step we consider the preconditioning matrix which is formed by the 
parameterized diagonal and lower triangular part of A 

R^V-f2C, (7) 

where 17 = diag{ujilm,uj 2 ln) with a;i,a ;2 G K. If i? is given by ([7]), then (jS]) 
becomes 


n{Ti,T2,uj2,a) = I - {V - QCy^TA 


or 

'H{Ti,T2,uJ2,a) = {V-nC)~\{I-T)V+{T - [2)C + TU] (8) 

and 

r]{Ti,T2,u>2,a) = {V - f2C)~^Tb. (9) 

Note that the parameter wi is absent in TL and rj. This is because the first m 
rows of C are zeros a fact which is carried over in matrix f2C also. 

The iterative scheme given by (0,®,® and ® will be referred to as the 
Generalized Modified Extrapolated SOR (GMESOR) method. In case a = 0 
this method was introduced in ® and proposed for further study. In the sequel 
to distinguish the dependence of GMESOR upon the parameter a we use the 
notation GMESOR(a). 

For {V — nC)~^ to exist we require 


det{V - n£) ^ 0. 

Because of 

Therefore, 

detiV — 17£) = (1 — aui 2 )^ det {A) det (Q) 0 


or 

auj2 ^ 1 ( 10 ) 

since the matrix A is symmetric positive definite and the matrix Q is nonsin¬ 
gular. In the sequel we require m to hold. 

The GMESOR(a) method has the following algorithmic form. 

The GMESOR(a) Method: Let Q € be a nonsingular and symmetric 

matrix. Given initial vectors and € K", and the parameters 

ti,T 2 yf 0, uj 2 ,a G with aui 2 yf 1. For k = 0,1,2,... until the iteration se¬ 
quence )^} is convergent, compute 


a;('=+i) = (1 - Ti)a;('=) -f riA-i(6i - 
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y{k+l) ^ y{k) ^ -Q-1 + (t 2 - W 2 )cc('=)] - T 262 I, 

where Q is an approximation of the Schur complement matrix A~^B. 

For special values of its parameters GMESOR(a) degenerates into known 
methods. Indeed, if w = ri = T 2 = a ;2 and o = 0 then GMESOR becomes 
the SOR-like method [H]; if w = ri = r 2 = a ;2 = 1 and a = 0 then it be¬ 
comes the preconditioned Uzawa method m-, and if ri = wi, r 2 = a ;2 and 
0 = 0, then it becomes the GSOR method [5]. By comparing the algorithmic 
structures of GMESOR(a) and GSOR, one can verify that the former has an 
additional matrix times a vector computation. Finally, if 

Ti = uj, — = 7 and — = T ( 11 ) 

1 — auj2 1 — auj2 

then the GMESOR(a) method becomes the Generalized Inexact Accelerated 
Overrelaxation (GIAOR) method [5] and if 

7"2 

Ti = uj,T 2 = UJ 2 and - = r ( 12 ) 

1 — aw 2 

the GMESOR(a) method becomes the Parametrized Inexact Uzawa (PIU) 
method [12] when P — A. The following theorem establishes the functional 
relationship between the eigenvalues A of the iteration matrix 'H{ti,T 2 ,uj 2 , a) 
and the eigenvalues /i of the associated matrix J = Q~^B^A~^B. 

Theorem 2.1 Let A £ symmetric positive definite, B £ 

of full column rank and Q £ be nonsingular and symmetric. If X ^ 1 —ti 

is an eigenvalue of the matrix'H{ ti,T 2 ,uj 2 , a) of the GMESOR(a) method and 
if p, satisfies 

+ + + = (13) 

\ 1 — aUJ2 J 1 — CL^2 

where auj 2 1, then p is an eigenvalue of the matrix J = Q~^B'^A~^B. 
Conversely, if p is an eigenvalue of J and if X ^ 1 — ti satisfies ilA) . then X 
is an eigenvalue o/'H(ti, r2, W2, a). In addition, A = 1 — ti is an eigenvalue of 
H{Ti,T 2 ,uJ 2 ,a) (if m > n) with the corresponding eigenveetor (x^,0)^, where 
X £ N'{B'^) and is the null space of B"'". 

Proof The eigenvalues p of the matrix J = Q~^B^ A~^B are real, positive and 
non-zero. Let A be a nonzero eigenvalue of the iteration matrix 'H(Ti,T 2 ,uj 2 ,a) 
defined in ([5|), and {x,y)'^ £ be the corresponding eigenvector. Then, 

'H(ri,r 2 ,a; 2 ,a) = A (14) 

or, from ([ 8 ]) we have 

[(j_T)p+(T_r2)£ + ™](^^^ =A(P-12£)(^^^. (15) 
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= A 


A 0 

—uj2B’^ (1 — au}2)Q 


From (I124|) and it follows that 
(1 - ti)A -tiB \ f X 

{t2-uj 2)B'^ {l-auj2)Qj\y 

Decoupling we obtain 

J (1 — ti)Ax — TiBy = XAx 

1 (t 2 — uj2)B'^x + (1 — auj2)Qy = —Xuj2B'^x + A(1 — aLj2)Qy 

or equivalently 

f (1 - Ti - X)x = TiA-^By 
I (1 - A)(l - auj2)y = [(1 - A)a ;2 - t2]Q~^B'^x. 

Multiplying the first equality in (11251) by Q~^B^, we obtain 

(1 — Ti — X)Q~^B^x = tiQ~^ B^ A~^ By, 

or, when X ^ 1 — ri, we have 


(16) 


Q-^B'^x = 


n 


-Q-^B^A-^By. 


1 — Tl — A 

From (I126p and the second equality in (11251) it follows that 

(1 - A)(l - aw 2 )(l - A - Ti)y = [(1 - A)a ;2 - T2]TiJy. 


(17) 


(18) 


If A = 1 — Tl 7 ^ 0, we have from (11251) that By = 0 and ti(1 — auj 2 )Qy = 
{tiuj 2 — T 2 )B'^x. It then follows that y = 0 and x S J\f{B'^). Hence, A = 1 — ti 
is an eigenvalue of H(ti , T 2 , W 2 , a) with the corresponding eigenvector {x’^ , 0 )^, 
where x G N[B'^). Therefore, because of (ITi^ . the eigenvalues A (except for 
A = 1 — Tl) of the iteration matrix 'H(ti, T 2 , W 2 , a) of the GMESOR(a) method 
and the eigenvalues y of the matrix J are related through the functional rela¬ 
tionship 

(1 - A)(l - aw 2 )(l - A - Tl) = [(1 - A)a ;2 - T 2 ]Ti/r, 

namely A satisfies the quadratic equation ((T5)) . □ 

From the above theorem we can obtain the following functional relationships 
for the GESOR(a), SOR-like(o) and GSOR(a) methods. 


Corollary 2.1 Under the hypothesis of Theorem, \2. 1\ 

1. The nonzero eigenvalues of the iteration matrix 'H{t, uj 2 , a) of the GESOR(a) 
method are given hy X = \ — t or if auj 2 ^ 1 by 


A^ + A T - 2 -f 


TUJ2 

1 — aw 2 ' 


y] +1-T- 


t(t - a;2) 

1 — aw 2 


^ = 0 . 


(19) 


2. The nonzero eigenvalues of the iteration matrix C{ui, a) of the SOR-like(a) 
method are given hy X = 1 — u: or if aui ^ 1 by 


X^ + Xluj-2 + 


1 — auj 


/ij “1-1 — LU — 0. 


( 20 ) 
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3. The nonzero eigenvalues of the iteration matrix C{iOi,uj 2 , a) of the GSOR(a) 
method are given hy \ = \ — uji or if auj 2 1 by 

+ x(uJi-2+ ^ ) + 1 - a;i = 0. (21) 

\ 1 - auj2 ) 

Proof The iteration matrix 'H(r, W 2 ,a) is obtained by letting r = ti = r 2 in 
TL(Ti,T 2 ,uJ 2 ,a). By following a similar approach as in the proof of Theorem 12.11 
we find the functional relationship dT51) . Similarly, we find (1^01) and □ 

Note that the above functional relationships are generalizations of the original 
SOR-like and GSOR methods. Indeed, if a = 0, then from (EHl) we obtain 
the functional relationship of the SOR-like method [21], whereas from (BID we 
obtain the functional relationship of the GSOR method [8]. □ 

Another preconditioning matrix R, which is formed by the upper triangular 
part of A is the following 

R = V- fiU. (22) 

Using (1^^ in ([5]) then ([5]) becomes the backward form of the GMESOR(a) 
method, which will be referred to as the Generalized Modified Extrapolated 
Backward SOR(a) (GMEBSOR(a)) method. From ([B]), because of (1^ . the 
iteration matrix of the GMEBSOR(a) method is given by 

/C(ri,T2,a;i,a;2,a) = I - {V - QU)~^TA 

or 

/C(Ti,T2,a;i,W2,a) = {V - - T)V + {T - f2)U + TC] (23) 

and 

fc(Ti, r 2 , wi, a; 2 , a) = (V — f2U)~^Th. (24) 

For {V — QU)~^ to exist we require 

det(T> - QU) ^ 0. (25) 

Because of 

Therefore, (ESI) becomes 

det()D — QU) = [1 — (1 — a)a; 2 ]" det A det Q 0 
or 

(1 - a)uj2 ^ 1 (27) 

since the matrix A is symmetric positive definite and the matrix Q is nonsin¬ 
gular. The GMEBSOR(a) method has the following algorithmic form. 

The GMEBSOR(a) Method: Let Q G be a nonsingular and sym¬ 

metric matrix. Given initial vectors € R'" and y^^^ € R", and the param¬ 
eters ti,T2 yf 0, wi, a; 2 , a G R with (1 — a)uj 2 yf 1. For fc = 0, 1 , 2, ... until the 
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iteration sequence 


y{k+l) _ y{k) _|_ 


T2 


1-(1 


^y{k) ^T| convergent, compute 
-& 2 ) 

- a)ijJ2 


a;('=+i) = (1 - {ri(6i - 


where Q is an approximate (preconditioning) matrix of the Schur complement 
matrix B^A~^B. 

As a by-product of the GMEBSOR(a) method we obtain the backward 
schemes corresponding to the GESOR(a) and GSOR(a) methods. For r = 
Ti = T 2 , we have the GEBSOR(a) method and for ti = wi and T 2 = W 2 we 
have the GBSOR(a) method. 

Corollary 2.2 Under the hypothesis of Theorem A‘2. 1\ 

1. The nonzero eigenvalues of the iteration matrix IC{Ti,T 2 ,uJi,u! 2 ,a) of the 
GMEBSOR(a) method are given hy \ = 1 — ti or if — a)uj 2 ^ 1 by 


X^ + X 



T2U}1 

1 — (1 — a)oj2 


T 


-I- 1 — Ti -I- 


T2{ti 

1-(1 


iXl) 

a)uj2 


y = 0. 


(28) 


2. The nonzero eigenvalues of the iteration matrix K(t, 101 , 1 x 2 , a) of the GEBSOR(a) 
method are given hy X = \ — t or if (1 — a)x 2 ^ 1 hy 


X^ + X 



TXl 

1 — (1 — a)x2 


T 


-I- 1 — r -I- 


r(r - xi) 

1 — (1 — a)x2 


^ = 0 . 


(29) 


3. The nonzero eigenvalues of the iteration matrix M.(xi, X 2 , a) of the GBSOR(a) 
method are given by X = 1 — xi or if (1 — a)x 2 ^ 1 by 


-|- A (xi — 2 -|- --—-- —-|- 1 — wi — 0. (30) 

V 1 - (1 - o)a;2 J 

Proof Following a similar approach as in the proof of Theorem 12.11 and using 
the iteration matrix IC{Ti,T 2 ,xi,X 2 ,a) given by (1^51) . we find the functional 
relationship (E51) . Similarly, we find the functional relationships (1^ and (15(11) . 

□ 

Note that the GMEBSOR(a) method has four parameters instead of three 
as the GMESOR(a) method whereas the GEBSOR(a) method has three pa¬ 
rameters instead of two as the GESOR(a) method. If a = 1, then (1501) becomes 
the functional relationship of the GSOR(a) method. 


2.2 Gonvergence 

In this section we develop the convergence analysis of the GSOR(a) and 
GMESOR method s as well as their corresponding backward counterparts. In 
particular, we derive sufficient conditions for the GSOR(a) and the GMESOR 
method to converge under the assumption that the eigenvalues of the matrix 
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J are all real. The sign of J’s eigenvalues depends upon the properties of the 
matrix Q. Specifically, we assume that the matrix Q is symmetric positive or 
negative definite. 

2.2.1 The GSOR(a) method 

The next theorem provides sufficient conditions for the GSOR(a) method to 
converge if the matrix Q is symmetric positive definite and a 7 ^ 0. 

Theorem 2.2 Let A G and Q G be symmetric positive definite 

and B G he of full column rank. Denote the minimum and the maximum 

eigenvalues of the matrix J = by Pmin and Pmax, respectively. 

Then p{L(u}i,uj 2 ,a)) < 1 if the parameters wi anduj 2 lie in any case of Table 

m 


Table 1 Sufficient conditions for the GSOR(a) method to converge if fimin > 0. 


Condition 

Cases 

CJ2 — Domain 

u)\ — Domain 

a > 0 

1 

2(2 — oji) 

0 < W2 < - ^ -T 

^if-i-max 2q,(2 ^ 1) 

0 < < 2 

a < 0 

2 


4a 

^^2 f \ \ 0 J 

^If-f^max 2(2(2 ^1) 

2a Id'Tnax 

3 

0 < UJ2 

4a 

0 < Wi < - 

2a Umax 

4 

2(2 — LOi) 

4a 

H“ 2(2(2 — UJl) 

2a l-^min 


Proof Recall lCorollarv l 2 . 1 |l that A = 1—wi 7 ^ 0 is an eigenvalue of £(a;i, a; 2 , a) 
and if A 7 ^ 1 — wi then the eigenvalues of C{uJi,uj 2 ,a) are given by (1^ . If 
A = 1 — wi 7 ^ 0, then the GSOR(a) method is convergent if and only if |A| < 1, 
that is |1 — wil < 1 , or 

0 < < 2. (31) 

If A 7 ^ 1 — wi and aw 2 1, then (1^ holds and by Lemma 2.1 page 171 of [5S] 
it follows that the GSOR(a) method is convergent if and only if 


where 


and 


|c| < 1 and | 6 | < 1 + c 


c = 1 — Wi 


b = 2 — LOi — 


UJ1OJ2P 


(32) 


(33) 


(34) 


1 — 0102 

From the first part of (11361) . because of (1551) . it follows that (1511) holds also in 
this case. From the second part of (11361) . because of (1551) and (1551) . it follows 
that 

UJiijJ2P 


2 — ijj\ — 


<2 — LVi 


1 — auj2 
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or 


0 < 


0 J 2 


< 


2 ( 2 -wi) 


(35) 


1 — 0102 Wi/i 

In order for to hold we distinguish two cases. Case I: 012 > 0 and 1 —0012 > 
0 and Case II: a ;2 < 0 and I — au )2 < 0. For each of theses cases we will 
distinguish two subcases, (i) a > 0 (ii) a < 0. In the sequel we will study 
the subcase (i) of Case I, since the other cases can be treated similarly. For 
subcase (i) of Case I 


1 


From (IMII . we have 


0 < a ;2 < 


+ 2n(2 — loi')^uj 2 ^ 2(2 — a^i). 


(36) 


(37) 


We distinguish two subcases: (ii) wi/i + 2a(2 — wi) >0 and (i 2 ) wi/i + 2a(2 — 
wi) < 0. In the sequel we will only treat subcase (ii) since the other case can 
be treated similarly. If ujifi + 2a(2 — wi) > 0 then 


4a > L0i{2a — /i). 

Next, we distinguish three subcases: (a) a > ^fimax (b) a < \^min (c) 
a < 2/1 max’ 

(a) o > h^max- From (l38ll we have 


(38) 


< 


Wi < 


4a 


2 a — fj,„ 

Combining OT and dSH), it follows that 


(39) 


0 < cji < min < 2, 


4a 


or 


Moreover, from (EZl we have 

0 < a ;2 < 

which, because of (1551) . becomes 


2 a l^mi 

0 < Wi < 2. 

2(2 -wi) 

^if^max 4” 2a(2 CUi) 


0 < 012 < min <( — 


2 ( 2 -wi) 


a LOl^J.rnax + ‘2a{2 - lOi) 


(40) 

(41) 

(42) 


which yields (IdTll again. Therefore, for case (a) we have that (l40ll and (1411) 
hold. 

(b) a < ^firnin- From (|551) we have 


UJi > 


4a 


2 a — firi 


(43) 
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Combining (l!m) and (HHll . it follows that 


max 


4a 


2a — fjLri 


< wi < 2, 


which yields (HUl) . Therefore, for case (b) we have that (pH)) and (HD) hold also 
as in case (a). 

(c) < a < Let Q;, /? be two positive integers such that = 

max{/x|/r < 2a}, = min{/x|/r > 2a}. Next, we distinguish two cases: (i) 

l^min — (ii) /^/3 ^ ^ ^max- 

Case (i): /imm < M < Ma- Following a similar approach as in Case (a), we have 
that (HOI) holds and 


0 < a ;2 < 


2 ( 2 -wi) 

4“ 2a(2 — cui) 


(44) 


Case (ii): fJ-p < fJ- < Umax- Following a similar approach as in Case (b), we 
have that (liO]) holds and 


0 < a ;2 < 


2 ( 2 -wi) 


^ll^max 4” 2a(2 CUi) 
Combining and (HHll it follows that 

2(2 — cui) 2(2 — cui) 


0 < 0^2 < min 


UJiHa + 2a(2 - Wi) ’ LOi^lmax + 2a(2 - Wi) J ’ 


(45) 


(46) 


which is equivalent to (HU- Hence, case 1 of table [T] is proved. Following a 
similar treatment we can prove the rest of the cases of Table [TJ □ 


Corollary 2.3 Under the hypothesis of Theorem \2.2\ a,nd if a = 0 then p(£(a;i, ^ 2 )) < 1 

if 

0 < wi < 2 and 0 < 102 < — -(47) 

^l^max 

Proof If we let a = 0 in (I^Tl) and follow a similar approach as in the proof of 
Theorem we can verify that (1471) holds. □ 

Note that (HTl) was also obtained in [5]. The following corollary gives sufficient 
conditions for the GBSOR(a) method to converge. 

Corollary 2.4 Under the hypothesis of Theorem \ 2. SA p{M{uJi,uj 2 ,a)) < 1 if 
the parameters wi and UJ 2 He in any case of Table [H 

Proof Using the functional relationship (IHHll and following a similar approach 
as in the proof of Theorem 12.21 we have 


0 < wi < 2 and 0 < 


UJ2 

1 — (1 — a)w2 


2 ( 2 -cui) 
ujip 


(48) 


Note that the second part of (H51) is the same as (1551) where now 1 — a appears 
instead of a. This occurs because the preconditioning matrix R is given by 
(1^^ and U is expressed in (|3]) in terms of 1 — a. Therefore, if we let 1 — o in 
place of a in Table [1] we obtain Table [H □ 
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Table 2 Sufficient conditions for the GBSOR(a) method to converge if firnin > 0. 


Gondition 

Cases 

UJ 2 — Domain 

cji — Domain 

a < 1 

1 

2(2 — oji) 

0 < a;2 < -^-^4--- 

H“ 2(1 u)(2 ) 

0 < oJi < 2 

1 < a 

2 

2(2-uii) 

4(1 — a) 

^ll-i'max 2(1 u)(2 ^l) 

4(1 Oj) l^max 

3 

0 < aj2 

4(1 — a) 

0 < wi < —-^- 

2(1 CL) f-L-max 

4 

2(2 — oji) 

4(1 - a) 

^1 P'TTiax 2(1 u)(2 

4( t Cl) l-l'min 


Corollary 2.5 Under the hypothesis of corollary \2.4\ and ifa = l then (wi, a; 2 ,1)) < 1 

if 

0 < LUi < 2 and 0 < uj 2 < — -(49) 

^ll^max 

Proof If we let a = 1 in (1481) then follows immediately. □ 

If the matrix Q is symmetric negative definite, then we have the following 
theorem. 

Theorem 2.3 Let A G be symmetric positive definite, B G be 

of full column rank and Q G be symmetric negative definite. Denote 

the minimum and the maximum eigenvalues of the matrix J = Q~^B^A~^B 
by Umin and pimax, respectively. Then p{C{uJi,u} 2 ,a)) < 1 if the parameters 
wi andu >2 He in the following cases of Table\^ 


Table 3 Sufficient conditions for the GSOR(a) method to converge if fimax < 0. 


Condition 

Cases 

(jj 2 — Domain 

u)\ — Domain 

a < 0 

1 

2(2 — oji) 

- ^^ -- < 4J2 < 0 

2fl(2 U-!l) 

0 < a;i < 2 

a > 0 

2 

~ ^ 2(2-uji) 

4a 

yO u ..J ^"^2 

H” 2(2(2 ) 

2a l-i'min 

3 

aJ2 < 0 

4a 

0 < tJi < -- 

2a l-^min 

4 

2(2 — oJi) 

4a 

H“ 2(2(2 U-!l) 

2a fJ'max 


Proof Using the functional relationship (I^Tl) and following a similar approach 
as in the proof of Theorem 12.21 taking into consideration that pmax < 0 we 
have 

0 < < 2 and < ——— < 0. (50) 

ojip 1 — au >2 

From (1501) the cases presented in Table |3] can be readily verified. □ 

Corollary 2.6 Under the hypothesis of Theorem \2.d\ a,nd if a = 0 then p{C{uJi,uj 2 )) < 1 

2(2 — a;i) 

0 < wi < 2 and - < ClI 2 < 0. 




( 51 ) 
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Proof Using the functional relationship m and following the proof of The¬ 
orem [5131 we have that if a = 0 in (IHHl) then (IHT1) follows. □ 

The above result was also obtained in [8]. 

Theorem 2.4 Let A G be symmetric positive definite, B G be 

of full column rank and Q G be symmetric negative definite. Denote 

the minimum and the maximum eigenvalues of the matrix J = A~^B 

by pL-min and Pmax, respectively. Then p{A4{uJi,uj2,a)) < 1 if the parameters 
uji andu >2 He in the following cases of Table\^ 


Table 4 Sufficient conditions for the GBSOR(a) method to converge if fimax < 0. 


Condition 

Cases 

ijJ 2 — Domain 

u)\ — Domain 

a < 1 

1 

2(2 — LOi) 

< U)2 <0 

^ll^min 2(1 <2')(2 

0 < oji < 2 

a > 1 

2 

~ ( 2(2 - 

4(1 — a) 

\ ^ / / \ / \ ^ 

H“ 2(1 u)(2 ^l) 

J- f-^min 

3 

UJ2 <0 

4(1 - a) 

2(1 n) f-^min 

4 

C 

\ 

c 

\ 

3 

1 

4(1 - a) 

^If^min H“2(l <2')(2 

■^( J- U) l^max 


Proof Using the functional relationship (1301) and following a similar approach 
as in the proof of Theorem 12.21 taking into consideration that Pmax < 0 we 
have 

0 < wi < 2 and —r— < 0. (52) 

uJip 1 — (1 — a)uj2 

From (1551) the cases presented in Table |3] can be readily verified. □ 

Corollary 2.7 Under the hypothesis of Theorem \2.4\ and if a = 1 then p(A4(uJi,uj2)) < 1 
*/ 

0 < wi < 2 and — -^ < a ;2 < 0. (53) 

Proof Using the functional relationship (1301) and following the proof of The¬ 
orem |2T4] we have that if a = 1 in (l52|) then (l53)) follows. □ 

In the sequel we study the convergence analysis of the GMESOR method under 
the same assumptions. 

2.2.2 The GMESOR method 

The next theorem provides sufficient conditions for the GMESOR method to 
converge if the matrix Q is symmetric positive definite and a = 0. The study of 
the case a ^ 0 follows a similar but cumbersome approach as it requires many 
cases to be examined. This study will not have any substantial contribution 
since the minimum value of the spectral radius of the GMESOR(a) method is 
independent of a fTheorem l2.1()l) . meaning that for, say a = 0, the GMESOR 
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method will attain the maximum rate of convergence. So, we are interested to 
find the convergence ranges of the parameters of the GMESOR(a) method for 
the simplified case when a = 0. 

Theorem 2.5 Consider the GMESOR method. Let A S and Q S 

be symmetric positive definite and B S ftg of full column rank. Denote 

the minimum and the maximum eigenvalues of the matrix J = A~^B 

by Umin and Umax, respectively. Then T 2 , W 2 )) < I if 


0 < Ti < 2, 0 < T2 < T2(^maa:) and ^2 < ^2 < ^2 ) (54) 


where 

4 1 _ 

l' 2 (/tmaa:) — , ^ 2 ^h'max') —'^2 and UJ 2 {t^max') 

'^if^max l^max 


2 - Tl T2 


'^ll^max 


2 ■ 

(55) 


Proof Recall that A = 1 — ti 0 is an eigenvalue of 'H{ti,T 2 ,uj 2 ) and if 
A ^ 1 — Tl then the eigenvalues of 'H(ti, T 2 , ^ 2 ) are given by (IT^ where a = 0. 
If A = 1 — Tl 0, then the GMESOR method is convergent if and only if 
|A| < 1, that is |1 — Ti| < 1, or 


0 < Tl < 2, 

(56) 

which is the first inequality of (IMl). If A 1— ti, then (TT^ holds and by Lemma 
2.1 page 171 of [SS], it follows that the GMESOR method is convergent if and 
only if (Il36|) holds where 

C = 1 - Tl +Ti(t 2 - W2)At 

(57) 

and 

6 = 2 — Tl — TlUJ2pi. 

From the first inequality of (Il36l) it follows that 

(58) 

0 < 1 + c < 2. 

(59) 


From the second inequality of (I136I1 . because of (EZD and (Ei), we have 

|1 + c — T 1 T 2 / 1 I < 1 + c 


or 


0 <^<l + c. 


Gombining (11391) and (1501) . it follows that 


0 < < l + c< 2. 


In order for m to hold we must have 


0<^<2 


(60) 


( 61 ) 


2 
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or, because of (IMll . 

0 < T2 < -, (62) 

Tl^ 

which proves the second inequality of (IMl) . Inequality (ICT) . because of (EH, 
becomes 

—-— <2 — Ti + TiT2fi — ria;2Ai < 2 
which is equivalent to 


1 

T 2 -< W 2 < 

A* 


2-Tl 

Ti/r 


II 

2 ■ 


(63) 


By studying the monotonicity of the right and left hand side of (1551) with 
respect to /i we obtain the third inequality of (IMll . □ 

The convergence conditions for GESOR are given by the following corollary. 

Corollary 2.8 Consider the GESOR method. Let A G and Q G 

be symmetric positive definite and B G R™^"- be of full column rank. Denote 
the minimum and the maximum eigenvalues of the matrix J = Q~^A~^B 
by Pmin and Umax, respectively. Then p{'H{t,uj 2 )) < I if 

0 < T < f{pmax) and U 2 (t) < UJ 2 < u> 2 (t), (64) 


where 


iAiir) = T- 



(65) 

( 66 ) 


Proof Letting r = ri = r 2 in (l54ll we obtain (l64l) . □ 

The convergence area for the GESOR method is illustrated in figure 1. Note 
that as Pmax increases the point of intersection of the two curves 0)2 (t) and 
W 2 (t) moves towards zero and the convergence area of the GESOR method 
shrinks. However, in practice Pmax usually is < 1. 

If the matrix Q is symmetric negative definite and a = 0 then we have the 
following theorem. 


Theorem 2.6 Consider the GMESOR method. Let A G R"^^™ be symmet¬ 
ric positive definite, B G R™^" be of full column rank and Q G R"^" be 
symmetric negative definite. Denote the minimum and the maximum eigen¬ 
values of the matrix J = Q~^B^A~^B by p-min and Pmax, respectively. Then 
P{U{ti,T 2,UJ2)) if 


0 < Ti < 2, T_ 2 {pmin) < ^2 < 0 and Ul 2 {pmin) < UJ 2 < 0 J 2 {Pmin)- (67) 


where 


X.2^L'min) — 


TlPr, 


^2^L'min) — 


, r 2 J _ , 

-h — and UJ 2 {Pn 


= r2- 


( 68 ) 
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X 


Fig. 1 Convergence area of the GESOR method for fimax = 0.99. 

Proof Following a similar approach as in the proof of Theorem l2.5l and using 
the functional relationship (1131) we can prove (1571) . □ 

Corollary 2.9 Consider the GESOR method. Let A € be symmetric 

positive definite, B G be of full column rank and Q G be symmetric 

negative definite. Denote the minimum and the maximum eigenvalues of the 
matrix J = by ytmin cindfimax, respectively. Then p{'H{t,uj 2 )) < 

1 */ 

(E.1 (.h"iTLin ) <C T 2, ^2 O'T^d fJ-min ^ 1; 

where 

/ \ ^ / 2 — T T j - f \ 1 

Z.lKh'min) — . ; — 1“ CLTlu, UJ2(/lmm) — ■ 

Y l^min l^min ^ h^min 

(70) 

Proof is proved by following a similar approach as in the proof of Theorem 
12.51 and using the functional relationship (ITOl) . □ 


2.3 Optimum parameters 

In this section we determine optimum values for the parameters of the iterative 
methods studied in the present section under the hypothesis that a 7 ^ 0 and 
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the eigenvalues of the matrix J are real. We assume that Q is a symmetric 
positive or negative definite matrix. 


2.3.1 The GSOR(a) method 


In the following theorem the optimum parameters for the GSOR(a) method 
are determined assuming that the matrix Q is symmetric positive definite. 

Theorem 2.7 Consider the GSOR(a) method. Let A £ and Q £ 

he symmetric positive definite and B £ he of full column rank. Denote 

the minimum and the maximum eigenvalues of the matrix J = 
by Pmin and pmax, respectively. Then the spectral radius of the GSOR(a) 
method, p(£(a;i, 012 , a)), is minimized for any a p-minV-max at 




iy/L min + Vp max)‘^ 

and its corresponding value is 


^ “ 1 “ -s/ 


\/M max - \fip n 
^/p max + ^/p n 


(71) 


(72) 


Proof The functional relationship (I^Tl) may be written as follows 


(A “b UJ\ — 1)('^ — 1) — — \uJ\Cj 2 p 


where 


W 2 

OJ2 = z - 

1 — auj2 


(74) 


with 0012 1. The optimum values of uji and 0)2 will be determined such 

that 

p(£(o;i, 0 J 2 , a)) = max |A| (75) 

lAmin 


is minimum. The real roots of (I73|) are the intersection points of the parabola 

(A + oil — 1)('^ ~ 1) 


9u,^ (A) = 


011012 


(76) 


and the straight lines 


^(A) — A/i, 0 < Pmin — P fi: Pvaax- (^^) 

Following a similar argument as in m page 111, h{X) are straight lines through 
the point (0,0) and guji{X) is a parabola passing through the point (1,0). The 
discriminant of (ED) is 

Z\(oil, 012 ,/!) = (2 - Oil - uJiU} 2 pf - 4(1 - oil). (78) 

Note that A{uji,u) 2 , p) < 0 for 0 < oii < dji{p) and A{oji,lj 2 , p) > 0 for 
dii{p) < wi < 2, where 

. , , ACj2P 
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If 0 < wi < wi(/^) then the value of p(£(wi,a; 2 ,«)) is 

|Ai| =|A^| = ( 80 ) 

where Ai and Xn are the two conjugate complex roots of (1^ . Furthermore, 
(IMl) is a decreasing function of wi. In case wi(^) < wi < 2 the roots of (1^ 
can be geometrically interpreted as the intersection of the curves gui^ (A) and 
h{X) = —XfJ,, as illustrated in figure 2, where we have assumed, without loss 
of generality, that h{X) = hi{X) = —Xg-max- The largest abscissa of the two 



points of intersection of h(X) and gujiiX) decreases with increasing uii. Indeed 
as uji increases, the intersection point (1 — a;i,0) of 5 cji(A) with the OX axis 
is moving towards to zero until guji{X) becomes tangent to h(X). Thus, for the 
fixed eigenvalue g of J, the value of wi which minimizes the zero of largest 
modulus of (1^ is wi(/r). Note that the straight lines /ii(A) = —Xfj,max and 
hwiX) = —Xgrnin include all the lines h{X) = —Xfi. Therefore, (17^ yields the 
two optima Cbi{gmax) and However, these values must be equal as 

there is only one optimum, hence 


max _ 4a)2A^min 
(1 + a}2/r max (1 4” 


(81) 


or 


UJ2 


1 


\/ gmingmax 


(82) 


which, because of CH), yields the optimum value for a ;2 given by the sec¬ 
ond part of CD). Substituting the value of W 2 in the expressions 0Ji{grnax) or 
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given by the first or second part of the equality (jST]), respectively, 
we obtain the optimum value of ui given by the first part of (ED- The spectral 
radius is given by 


p{C[uJi,uj 2 ,a)) = max{|Ai|, IAatI} (83) 

where Ai, Aat are the abscissas of the points of tangent of hi (A), hAr(A), respec¬ 
tively. For the minimization of p(£(a;i, a; 2 , o)) with respect to wi we require 

|Ai| = IAatI 


or 

Al = —Xn = (1 — , (84) 

where the last equality holds by the fact that Ai, Xn are the abscissas of the 
tangents hi (A) and hN{X), respectively. From (IMl) and (IMll it follows that 

p{£{uji,uj2,a)) = (1 

which, because of (ITT]) , yields ([72]). □ 

Theorem 12.71 finds the optimum values of the relaxation parameters wi and u }2 
of the GSOR(a) method. Letting a = 0 in dTTll we obtain the optima found 
also in [S]. Note that the parameter a has no impact on the spectral radius 
of the GSOR(a) method as one might have expected. The algebraic approach 
in [8] is similar to the one followed by |55] for determining the optimnm of 
the sole parameter in the SOR method. In case of GSOR(a), which has two 
parameters, there is an alternative less tedious algebraic approach (see [SS] pp. 
279-281). However, it remains to be verified whether either approach can be 
used to solve the problem of determining the optimum values of more than two 
parameters as is the case for the GMESOR(a) method. Our approach follows 
the geometric approach of Varga m for the determination of the optimum 
value of the parameter w in SOR. It should be noted that this approach is also 
mentioned in m but without a proof. 

Corollary 2.10 Consider the GBSOR(a) method. Under the hypothesis of 
Theorem \2. 7| the spectral radius of the GBSOR(a) method, p{M.{uji,u} 2 ,a)), is 
minimized for any a ^ 1 + y/PminPmax at 


a^lopt = 


4 .^/ P'minPn 


(v^ min + a/M max)^ 


and uj2„„t = 


(1 a) ^ P'minPn 


(85) 


and its corresponding value is 




max - Vu min 
y/U max + y/U min 


(86) 
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Proof We remark that the functional relationship (l!TO of GBSOR(a) is the 
same as that of the GSOR(a) method (I^T|) with the only difference that now 
we have 1 — a instead of a. Therefore, we have the same results as in Theorem 
1^:71 if we simply replace a with 1 — a. □ 

If the matrix Q is symmetric negative definite, the optimum parameters and 
the minimum spectral radius for the GSOR(a) method are given by the fol¬ 
lowing theorem. 

Theorem 2.8 Consider the GSOR(a) method. Let A G symmetric 

positive definite, B G ]R™xn of full column rank and Q G be symmetric 

negative definite. Denote the minimum and the maximum eigenvalues of the 
matrix J = Q~^A~^B by pimin and Pmax, respectively. Then the spectral 
radius of the GSOR(a) method p(£(a;i, W 2 ,«)), when the matrix J has negative 
eigenvalues, is minimized for any a y/Uminh-max at 


a^iopt = 


d-y/ P'minh'max 


(vIm min I + max I 


UJ2o^t = 


^ \/ l^minf^max 


and its corresponding value is 


(87) 




\/\t min I - \/\t ri 

'At min I + \/\t n 


( 88 ) 


Proof In this case /r < 0. Following a similar approach as in Theorem 12.71 
we obtain dUl) and (1881) . □ 

Under the hypothesis of Theorem 12.81 and if a = 0, these results were also 
obtained in [8]. 

Theorem 2.9 Consider the GBSOR(a) method. Let A G ]R™xm symmetric 
positive definite, B G R™xri ^g qJ column rank and Q G be symmet¬ 

ric negative definite. Denote the minimum and the maximum eigenvalues of 
the matrix J = Q~^B^A~^B by p,min and Umax, respectively. Then the spec¬ 
tral radius of the GBSOR(a) method p(M.{u}i,u) 2 ,a)), when the matrix J has 
negative eigenvalues, is minimized for any a ^ 1 — y/PminPmax at 


a^iopt = 


^y/ pminpmax 


{At min I + vIm max I 


and = 


1 a yj PminPmax 


(89) 


and its corresponding value is 


p{C{uJi 

opt 5 aJ2opt,a)) 


\/\t min I - '/\t ri 

\/\t min I + 'At Ti 


(90) 


Proof We remark that the functional relationship (15(11) of GBSOR(a) is the 
same as that of the GSOR(a) method (I5T|) with the only difference that now 
we have 1 — a instead of a. Therefore, we can apply the results of Theorem l2.8l 
by replacing a with 1 — a. □ 
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2.3.2 The GMESOR(a) method 


In the sequel we determine the optimum parameters for the GMESOR(a) 
method. 

Theorem 2.10 Consider the GMESOR(a) method. Let A G and Q G 

be symmetric positive definite and B G be of full column rank. De¬ 

note the minimum and the maximum eigenvalues of the matrix J = 
by tLmin and pLjnax, respectively. Then the spectral radius of the GMESOR(a) 
method, p{'H(Ti,T 2 ,uj 2 ,a)), is minimized for any a ^ —y/PminPmax at 


— T 2 „ 


PminPn 


( fT, -— rr, - ')2 

\-\/ f-^min I y H'max J 

and its corresponding value is 


and = = 

^ “T Y f^minl-^max 




y/r max - y/r D 
max + Vt n 


(91) 

(92) 

(93) 


Proof The functional relationship of the GMESOR(a) method is given by 
(fT^ or 

+ = (M) 

1 — aLU2 

The optimum values of ti , T 2 and UJ 2 will be determined such that 


piH{Ti,T2,uj2,a)) = max |A| 


(95) 


is minimum. The real roots of (USD are the intersection points of the parabola 

(A + ri-l)(A-l) 


gW = 


n 


and the straight lines 


_ “^2 - T2 - Aw2 

h(A) — /i, 0 <C fJ^rnin ^ ^ Mn 


(96) 


(97) 


1 — aui2 

Following a similar argument as in m page 111, h{X) are straight lines through 

( UJ2 — T2 \ 

0 ,- p and 5 ti(A) is a parabola passing through the points 

1 — aLJ2 J 

(1,0) and (1 — ri,0) (see figure 3). The spectral radius is given by 

p{'H{Ti,T 2 ,uj 2 ,a)) = max{|Ai|, |Ajv|} (98) 

where Ai,AAr are the abscissas of the points of tangent of /ii(A),/lAr(A), re¬ 
spectively, where now hi{X) = {ui2 — T 2 — Xui2)pmax and h^iX) = (a ;2 — T 2 — 
Xuj 2 )pmin- Therefore, 


|Ai| = ( 1 - n - —Pmax 

1 — auj2 


1/2 


(99) 





















Fig. 3 Graphs of (A), hi (A) and h]\i{\) in case the roots of (1941 1 are real. 


and 

IAatI = - Ti + ■ (100) 

\ i — au)2 J 

From (|98D it follows that the minimum value of r 2 , a; 2 , a)) is attained 

when 

|h| = lAjvl (101) 

which, because of (IMll and (IIOOII . implies 

UJ2= T2. ( 102 ) 

In case Ai and \n are the two conjugate complex roots of (|M)) . it follows that 
(I101|l must also hold for p{'H{Ti,T 2 ,u} 2 ,a)) to be minimized. So, (11021) holds if 
either (|94|) has real or conjugate complex roots. However, if (11021) holds, then 
m becomes 

A^ + A (ri - 2 + tiT2p) + 1 - n = 0, 


which is the functional relationship of the GSOR with 
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Therefore the optimum values of ri and T 2 are given by Wopt and Topt of 
[5], respectively, whereas the minimum value of p('H(ri, T 2 , W 2 , a)) is given by 
p{T-L{ujopt, Topt)) of [ 8 ]. Finally, using (11031) we find (l9^ . □ 

So, for the optimum values of its parameters, GMESOR(a) degenerates to the 
GSOR(a) method. 

Corollary 2.11 Consider the GESOR(a) method. Let A G q g 

be symmetric positive definite and B G of full column rank. De¬ 

note the minimum and the maximum eigenvalues of the matrix J = 
by Pmin and Pmax, respectively. Then the spectral radius of the GESOR(a) 
method, p('H{T,uj 2 ,a)), is minimized at 


^2opt ~ Topt 


and 


_ 4^ minl^max 1 - - 

'^0‘pt — ~ j ^ ^ j ! ^ ^ Y 2 ’ ^opt — yj j^maxl^r 

'^opt 


(Vt min + y/fi max)^ 

and its corresponding value is 

p{'H{Topt, aj2ppf, Oopt)) = 


max - y/r min 
y/r max + y/fi min 


(104) 

(105) 


(106) 


Proof Recall that GESOR(a) is obtained by setting ti = T 2 in GMESOR(a). 
Therefore, (nnn) and (I105|) are obtained by (IMI) and dMl), respectively, where 
now we require Topt = = T 2 „pt. □ 


Corollary 2.12 Consider the GMEBSOR(a) method. Under the hypothesis of 
Theorem \2. 1 ()\ the spectral radius of the GMEBSOR(a) method, p{K{Ti,T 2 ,uJi,uj 2 ,a)), 
is minimized at 


^lopt Tlppf , 


(107) 


where 


^lout 




(y/r min + y/r max 


and its corresponding value is 


and T 2 „^, 


1 — (1 — a)u)2 

y/ PminPmax 


(108) 


p(/C(ti 

opt 1 T2apt,^lopt,^2, a)) — 


y/fi max - y/fi D 
y/fi max + y/fi D 


(109) 


Proof Following a similar approach as in Theorem 12.101 using the functional 
relationship (l28ll and requiring |Ai| = IAatI we find 


UJl=Tl. (110) 

Therefore, (1^ because of (IllOl) becomes 

A^ + A (ti — 2 + TlT2/i) + 1 — Tl = 0 (111) 


T2 

1 — (1 — a)uj2 ’ 


with 


T2 = 


( 112 ) 
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which is the functional relationship of the GSOR(a) method (see (EH)) with 
the only difference that now we have 1 — a instead of a in f 2 , hence (110811 and 
(I109|l hold because of Theorem 12.71 □ 

Note that although the GMEBSOR(a) method has four parameters instead of 
three as in the GMESOR(a) method, both methods have the same minimum 
spectral radius. 


Corollary 2.13 Consider the GEBSOR(a) method. Under the hypothesis of 
Theorem, \2.10\ the spectral radius of the GEBSOR(a) method, p{IC{T,uj 2 ,a)), 
is minimized at 

^lopt ~ '^opt (113) 

and 


'^opt — 




min + max)^ 


— 


1 '^opt -sj prainPmax 
1 — 0 


and its corresponding value is 


p(/C(Topt, , o)) — 


y/U max - \/m min 
y/U max + y/J^ min 


(114) 


(115) 


Proof Recall that GEBSOR(a) is obtained by setting ri = T 2 in GMEBSOR(a). 
Therefore, (11131) and (11141) are obtained by (11071) and (110811 , respectively, where 
now we require Topt = = T 2 „^, . □ 

Our analysis so far shows that all the studied iterative methods (GMESOR(a), 
GMEBSOR(a)) have also the same rate of convergence as the PCG method for 
the optimum values of their parameters (see Theorems 12.7112.101 and corollary 

[HH). 


3 The Generalized Modified Preconditioned Simultaneons 
Displacement (GMPSD) method 

The Preconditioned Simultaneous Displacement (PSD) method was intro¬ 
duced in m- When the coefficient matrix A is two-cyclic the Modified PSD 
(MPSD) method was studied in [53], [55] ■ Motivated by our previous work 
we introduce the Generalized Modified PSD (GMPSD) method and study its 
convergence rate for the numerical solution of the augmented linear system 

O-ED- 


3.1 The functional relationship 

In the sequel, we let the preconditioning matrix R be the product of the 
lower triangular part with the upper triangular part of A in an attempt to 
obtain a better approximation of A and consequently an increase in the rate 
of convergence of the corresponding iterative method. Let 

R={V- nC)V-^{V - QU). 


(116) 





























26 


M. A. Louka et al. 


From (|ni) and (111611 it follows that the iteration matrix of ([S]) now is 

g{Ti,T2,uji,uj2,a) = I - {V - nU)-^D{V - f2£)-^TA (117) 

whereas t 7 (ti,T 2 ) in (jS]) corresponds to 

7(Ti,T2,a;i,W2,a) = {V - D{V - SlCy^Tb. (118) 


Note that this method has four parameters ti , r 2 , wi and uj 2 instead of three 
in the GMESOR method. The iterative scheme given by (O, (I117p and (111811 
will be referred to as the Generalized Modified Preconditioned Simultaneous 
Displacement (GMPSD) method. For {V — fiU)~^D{'D — J7£)~^ to exist we 
require 

det[{V-QC)V'^'^{V-nU)]^0. (119) 

Because of O 


R = {v-nc)v^^{v-m) 


Therefore, 


( ^ 

\—uj2B'^ (1 — aa;2)[l — (1 


ojiB 

— a)ui2]Q—uJiUJ2B'^ A~^B 

( 120 ) 


det(T' - QC)V~^{V - nU) = (1 - auj 2 T[l - (1 - a)a; 2 ]" det {A) det (Q) ^ 0 


or 

a 7 ^ - and uj 2^2 (121) 

since the matrix A is symmetric positive definite and the matrix Q is nonsin¬ 
gular. The GMPSD method has the following algorithmic form. 

The GMPSD Method: Let Q G be a nonsingular and symmetric ma¬ 

trix. Given initial vectors € M™ and € R", and relaxation faetors 
D, T 2 0, wi, a; 2 , a G R with a ^ and uj 2 yf 2. For fc = 0,1, 2,... until the 

iteration sequence )^} is convergent, compute 

{Bn{T2 - ria;2)a:W + T,u:2A-\b, - By«)]-r 262 } 

= (l-Ti)xW A-i {B [(wi - Ti)yW - Wiy('=+1)] + Tibi], 

where Q is an approximation of the Schur complement matrix B"^A~^B. 

Note that in the above algorithm we first compute and then 

whereas in the GMESOR method we had the reverse computations. If r = 

Ti = T 2 and w = wi = a ;2 we have the GPSD method. 

If a ;2 = 0 then the algorithmic form of the GMPSD method simplifies to 

y{k+i) = y{k) + r2Q-\B^x^’^^-b2) 

^{k+i) = +riA-^{bi-By^’^+^l) ^ ’ 

The above form is the same as that of the GSOR method if we use V — ilU 
instead of 27— fiC as the preconditioned matrix in the GSOR method and will 
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be referred as the simplified GMPSD method. In the following theorem we find 
the functional relationship for the GMPSD method between the eigenvalues A 
of the iteration matrix C/(ri,T2,a;i,W2,a) and the eigenvalues /r of the matrix 
J. 

Theorem 3.1 Let A G symmetric positive definite, B G 

of full column rank and Q G be nonsingular and symmetric. If X 7^ 1 —ti 

is an eigenvalue of the matrix Q{ti,T 2, to 1,102, a) and if pL satisfies 


A + A ( Ti — 2 + 


ria;2 + T2L0i — T 1 LJ 1 UI 2 
(1 — aaj2)[l — (1 — a)oJ2]' 


At +1 - T1 + 


n'r 2 — ti(jJ 2 — T2IV1 + T1UJ1UJ2 

(1 — aoJ 2 )[^ — (1 — a,)u) 2 \ 


At = 0, 


( 123 ) 

where a ^ \ and LO2 7^ 2 , then fi is an eigenvalue of the matrix J = Q~^A~^B 
Conversely, if p, is an eigenvalue of J and if X ^ 1 — ti satisfies \ 12 S\) . then X 
is an eigenvalue of 0 {ti,T 2 ,uji, 102,0). In addition, X = 1 —ri is an eigenvalue 
of G{ri,T2,L0i,L02,a) (if m > n) with the corresponding eigenvector {x'^, 0 )'^, 
where x G Af{B'^). 


Proof Glearly, the eigenvalues fj, of the matrix J = Q~^B^A~^B are real and 
non-zero. Let A be a nonzero eigenvalue of the iteration matrix G{ti, T2,loi,uj 2, a) 
defined in ( 1113 , and {x, G be the corresponding eigenvector. Then, 

we have that 

g{Ti,T 2 ,U!i,L 02 ,a) 

or because of (Uni) 



[{V - nC)V-^{V - QU) - TA] 



x{v-nc)v-\v-nu) 



From (I 124 L because of o, we have that 


{1 -ti)A {ooi-ti)B 

(t2 — oj2)B'^ (1 — aw2)[l — (1 — a)oj2]Q — ooiui2B'^A~^BJ \y 

— X f ^ oJiB 

—uj 2B'^ (1 — 0012 )[1 — (1 — a)(j02]Q — LO1LO2B"’"A~^B 

Decoupling we have that 



{ (1 — ti)Ax + (uji — Ti)By = XAx -\- XuoiBy 
(t2 — uj2)B'^x + {(1 — aa;2)[l — (1 — a)u}2]Q — ujiuj 2B'^ A~^B}y 

= —Xlo2B'^x + A {(1 — ai02)[l — (1 — a)i02]Q — ujiL 02B'^ A~^ B}y 

or equivalently 


(1 — Ti — A)a; = [(A — l)wi -I- ti\A ^By 

(t2 — L02 + Xlo2)Q~^B^x = (A — 1 ) {(1 — aw2)[l — (1 — 0)^2]/ — a;iu;2^} y- 

( 125 ) 
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From the first equality in (112511 we get 

(1 - n - X)Q~^B^x = [(A - l)a;i + n] Jy, 
and hence, when \ ^ 1 — ti, 

Q~^B'^x = —— ^ jy ^226) 

1 — Ti — A 

It then follows from (11261) and the second equality in (I125p that 

(A - 1)(1 - aw2)[l - (1 - a)w2](l - n - A)y 

= {[(A - l)w2 + ''■ 2 ][(A - l)a;i + n] + (A - 1)(1 - n - A)a;ia;2}Jy. 

If A = 1 — Ti 0, then from the first and the second equality of (11251) we have, 
respectively. By = 0 and ri{(l — aw 2 )[l — (1 — a)u} 2 ]Q — ojilo 2 B'^A~^B}y = 
{tiuj 2 — T 2 )B'^x. It then follows that y = 0 and x G M{B'^), where M{B'^) 
is the null space of the matrix Hence, A = 1 — ti is an eigenvalue of 
t/(Ti,r 2 ,wi,072,a) with the corresponding eigenvector (a;^,0)^, where x G 
Af{B^). Therefore, the eigenvalues A (except for A = 1 — ri) of the matrix 
t/(Ti, r 2 , wi, 072, a) and the eigenvalues y of the matrix J satisfy the functional 
relationship 

(A — 1)(1 — auj2)[l — (1 — a)w2](l — n — A) 

= {[(A - 1)W2 + T2][(A - l)wi + Ti] + (A - 1)(1 - Ti - \)l0iL02} H- 
This means that A satisfies the quadratic equation (I123I1 . □ 

Corollary 3.1 Let A G be symmetric positive definite, B G be 

of full column rank and Q G R"^” be nonsingular and symmetric. 

1. The nonzero eigenvalues of the iteration matrix Q(t,uji,uj 2 , a) ofthe GMPSD(3) 

method are given by \ = 1 — t or if a ^ and 102 2 by 

\2 \ f n ™ \ 1 t(t — uj) 

^ (1 - auj 2 ){l - (1 - a)uj 2 ] V (1 - auj 2 )[l - (1 - 0 )^ 2 ] ^ ^ ^ 

(127) 

where 

Cj = UJi + UJ2 — (128) 

2. The nonzero eigenvalues ofthe iteration matrix S { 0 ) 1 , 102 , a) ofthe GMSSOR 
method are given by X = 1 — io or if a ^ and 0 J 2 by 

where uj is given by IMS)- 

The nonzero eigenvalues ofthe iteration matrix Q{t, 10 , a) ofthe GPSD method 
are given by X = 1 — t or if a ^ and lo 2 by 

o / TOO \ t(t — w) 

A +A ( r - 2 + —-g---- —y | +l-r+—-----= 0 

\ (1 — aa;)[l — (1 — a)a;J J (1 — oa;)[l — (1 — a)a;J 

( 130 ) 

















A comparison of the ESOR and PSD methods for augmented linear systems 


29 


where now 

u) = a ;(2 — (jj). ( 131 ) 

and 

aoj ^ 1 and (1 — a)uj ^ 1 (132) 

4- The nonzero eigenvalues of the iteration matrix S{uj,a) of the GSSOR 
method are given by X = 1 — vj or if a ^ and uj 2 by 

+ A — 2 + —-r-—— -+ 1 — w = 0 (133) 

\ (1 — aw)[l — (1 — ajwj / 

where uj is given by 

Proof The iteration matrix C/(r,a;i,a; 2 , a) of the GMPSD(3) is obtained 
by letting r = ti = r 2 in C/(ti, r 2 , wi, W 2 ,a) given by (111711 . Using the matrix 
t/(Ti, T 2 , wi, W 2 , a) and following a similar approach as in the proof of Theorem 
13.11 we find the functional relationship (11271) . Similarly, we find (I129p . (11301) 
and (|133l) . □ 


3.2 Convergence 


If the matrix Q is positive definite and a = 0 sufficient conditions for the 
GMPSD method to converge are given by the following theorem. 


Theorem 3.2 Consider the GMPSD method. Let A G and Q G 

be symmetric positive definite and B G be of full column rank. Denote 

the minimum and the maximum eigenvalues of the matrix J = Q~^A~^B 
by Umin andyimax, respectively. Then, p{Q{ti,T2,uji,uj2)) < 1 if the parameters 
Ti,T 2 ,uJi anduj 2 He in the region defined in the cases of Table\^with 0 < ri < 2 
and 




ti{2uj2 - T2) (ti - 2)(1 - UJ2) 1 

2(tiu;2 —T 2 ) T 1 UJ 2 — T 2 p 

ti{uj 2 - T2) Ti(l - a;2) 1 

Tia;2 — T2 T 1 UJ 2 — T2 p' 


-^21 


Tl 

5 

n 


W22(m) = 1 - 


T 1 T 2 PL 

4 


(134) 


Table 5 Sufficient conditions for the GMPSD method to converge. 


Cases 

UJ 2 — Domain 

uji — Domain 

T 2 — Domain 

1 

6^21 ^2 ^22 

^11 ^1 

4ri 

0 < ^2 < 2 

4 H“ T-| f-tmax 

2 

Ul2 < 


3 

L02 < 

4ri 

. , 9 < '^2 

4 + T., (Trnin 

4 

1 < LD2 ^22 

T2 < 0 
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Proof Recall that A = 1 — ri 7^ 0 is an eigenvalue of G{ti,T 2 ,uji,u} 2 ) and 
if A 7^ 1 — n then the eigenvalues of G{ti,T 2 ,uji,uj 2 ) are given by (|12dD . If 
A = 1 — Ti 7^ 0, then the GMPSD method is convergent if and only if |A| < 1, 
that is 11 — Ti I < 1, or 

0 < n < 2. (135) 

If A 7^ 1 — Ti, then (I123|) holds and by Lemma 2.1 page 171 of [55], it follows 
that the GMPSD method is convergent if and only if 

|c| < 1 and |6| < 1 + c (136) 


where 


and 


c = 1 — Tl 


ria;ia;2 — ria;2 — T2UJ1 + T1T2 
I — UJ2 


-M 


0 = 1 + c — -- 

1 — ClI 2 

From the first inequality of (11361) it follows that 

0 < 1 + c < 2. 


(137) 

(138) 


(139) 


From the second inequality of (I136I1 . because of (I138|) . we have 




2(1-a;2) 

Combining (I139p and (|140l) it follows that 


0 < , < 1 + c < 2. 


2(1-W2) 

In order for (11411) to hold we must have that 


or because of (11351) 


2(1 — W2) 


T2 4 

0 < —< 


1 — a;2 Tifi 

Inequalities (11411) . because of (11371) . become 


(140) 


(141) 


(142) 


ri(2a;2 - T2)fi 

2(1-W2) 


+ Tl — 2 < Wl 


T1W2 — T2 

1 — a;2 


At < Tl + 


Tl(aJ2 - T2)fJ, 
1 — UJ2 


(143) 


In the sequel we distinguish the following two cases to study Case I: 

T2 > 0 and 1 — W2 >0 and Case II: T2 < 0 and 1 — 012 < 0. In addition, 
we distinguish the following two subcases for each of the above cases, (i): 
Tia;2 — T2 > 0 and (ii): T1W2 — T2 < 0. Next, we will study only the subcase (i) 
of Case I, since the other cases can be treated similarly. For this case, we have 
that 
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and from the second part of p42D 


a;2 < 1 


TlT2^J, 

4 


(145) 


From (I144p and (11451) it follows that 


— < U }2 < min 11,1-1, 0 < T 2 < Ti 

Tl 14 1 


or 


T2 

— <UJ2<1 
Tl 


TlT2fJ. 

4 


0 <T2 < Tl 


(146) 


which holds if ^ < 1 — Therefore, we have that (11461) holds if 


4ti 

^21 <UJ2 < ^ 22 ( 1 ^), 0 < r2 < 2 (147) 

4 + Tj /r 

where a;2i,W22(M) given by (I134|) . Furthermore, from (I143|) . we have that 

a;*i(/r) < < Wi2(Ai) (148) 

where a;*]^(/x), a;j['2(/^) are given by (11341) . Studying themonotonicity ofw22(At)j<^ii(M) 
and uJi 2 ifJ') with respect to ^ we have that sign = —1, sign = +1 

and sign = +1- Hence, case 1 of Table is proved. Treating similarly 

subcase (ii) of Case 1 and subcases (i) and (ii) of Case 11, we can prove the 
rest of the cases in Table [5] □ 

The convergence conditions for the GMPSD(3) are given by the following. 


Corollary 3.2 Consider the GMPSD(3) method. Let A € and Q G 

l^nxn symmetric positive definite and B G be of full column rank. De¬ 

note the minimum and the maximum eigenvalues of the matrix J = Q~^A~^B 
by pLmin and g.max, respectively. Then, p(C/(r,a;i, W2)) < 1 */ 

0 < T <2, LV2 ^ aS2^LTn.ax) and UJ^^(^g.rnax) ^ OJi < (/^max) (149) 


where 


((A^) = l- 


^13 


(m) 


T — U)2 
l — i02 




2—t _j_ r — 2(jj2 
rp 2(1—u;2) ’ 


(150) 


Proof Letting a = 0 in the functional relationship (11271) and following a similar 
approach as in the proof of Theorem 13.21 we can prove (11491) . □ 

Note that analogous results hold when Q G is symmetric negative 

definite. 
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3.3 Optimum parameters 


In the following theorem the optimum parameters of the GMPSD method are 
determined assuming that the matrix Q is symmetric positive definite and 
o ^ 0. 


Theorem 3.3 Consider the GMPSD method. Let A G and Q G 

be symmetric positive definite and B G R™^"- be of full column rank. Denote 
the minimum and the maximum eigenvalues of the matrix J = Q~^A~^B 
by Umin and g^max, respectively. Then the spectral radius of the GMPSD method, 
p{Q{Ti,T 2 ,uii,uj 2 ,a)), is minimized for any uj 2 at 


OJlcnt - 


IXmin Pmax 

“ ( / ,, . -L / ,< ^2 

\yj h^min \ y h^max ) 

and its corresponding value is 


n 

opt opt - aJ2) 

(1 — aw2)[l — (1 — a)u}2] 


and T 2 „„, = 




Pminpmax 


y/P max - y/fi min 
y/P max + y/P min 


(151) 

(152) 


(153) 


Proof Following a similar approach as in Theorem 12.101 using the functional 
relationship (11231) and requiring |Ai| = |Aw| we find 


Tiwia;2 — r2a;i — ria;2 + tiT2 = 0. (154) 

Therefore, (112311 because of (I154|) . becomes 

A^ + A (ri — 2 + Tit2p) + 1 — Ti = 0 (155) 


with 


T2 


_ T2 _ 

(1 — aa;2)[l — (1 — 0)^2] 


(156) 


which is the functional relationship of the GSOR method [5] with the only 
difference that now we have (1 — 0022) [1 ~ (1 ~ 0)^2] instead of 1 — aui 2 in the 
denominator of T2 (see (| 1031) 1. hence (I151L follows from (115411 whereas (115211 
and (115311 hold because of (I155L (11561) and Theorem 4.1 in [5]. □ 


Corollary 3.3 Consider the simplified GMPSD method. Let A G R™^™ and 
Q G R"^" be symmetric positive definite and B G be of full column 

rank. Denote the minimum and the maximum eigenvalues of the matrix J = 
Q~^B^A~^B by Pmin and Pmax, respectively. Then the spectral radius of the 
simplified GMPSD method, p(1/(ti, T 2, wi, 0, 0)), is minimized at 


^lopt Dopf> 


(157) 


l^minl-^max 


{y/P min + y/P n 


and T2„p, = — 
yfiP 


minp^max 


( 158 ) 
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and its corresponding value is 


opt 5 T2 

Opt 5 Wi 

opt 5 0,0)) = 


\/M max - yjj- min 




^/Pr, 


(159) 


Proof Letting uj 2 = 0, (I151|) . (11521) and (11531) yield (jl57l) . (I158p and (I159p . 
respectively. □ 

It is worth noting here that the optimum values of and T 2 ^^^ of the sim¬ 
plified GMPSD method are identical to the optimum values of and u}2apt 
of the GSOR method, respectively. 


Corollary 3.4 Consider the GMPSD(3) method. Let A G ]^”rxm q ^ 
be symmetric positive definite and B G ftg of full eolumn rank. De¬ 

note the minimum and the maximum eigenvalues of the matrix J = 
by Pmin and Pmax, respectively. If Pmax < \ or if Pmax > \ and either (i) 
Pmin < P* or (a) prriin > P* and Oi < tt < 02, then the spectral radius of the 
GMPSD(3) method, p{Q{T,uJi,uj 2 ,a)), is minimized at 


Topt 0J2p t 

a^lopt = -j-: 

1 - 0J2ppt 


( / 7 ^+ /jT- V = 

\\ h^min I s/ H'max J 


and its eorresponding value is 

p{g{Topt,UJi^p,,UJ2ppt,a)) = 

where 


P = 


l-^r] 


(1-20W)2’ 


Oi = 


1 ± \/I — a(l — a)a 

yfP max - Vp min 
y/P max + y/P min 

2 

02 = 


with 


a = 4(1 — M) and M = 


a—\Ja[p — 4) cr+y/ a(a — 4) 

^Pminpn 


{y/P min + yfp maa;)^ 


(160) 

(161) 

(162) 

(163) 

(164) 


Proof Recall that GMPSD(3) is obtained by setting ti = T 2 in GMPSD. 
Therefore, (UlOl), mn) and (11621) are obtained by (11511) . (11521) and nisi) re¬ 
spectively. In particular, by letting = T 2 ppt it follows from (11521) that 


o(l — a)uj 2 — 0 J 2 + I — M = Q, 

where M is given by (Ell). This quadratic has real roots when 

a^a — aa 1 > 0, 


(165) 

(166) 


where cr is given by (11761) . Gonsidering (11781) as a quadratic we distinguish two 
cases. Gase 1: Z\a < 0, Case 2: Z\a > 0 where Aa = cr{a — 4). 




















































34 


M. A. Louka et al. 


Case 1: Aa < 0. In this case we require cr > 0 since cr — 4 < 0 or in view of 

(CZSD 

min (1 - 2^ 

max ') > -\/m max ■ (167) 

But, (I179P holds if either fj^max < -j or if ^j^max > j and /imm < fJ-* and (i) is 
proved. 

Case 2: Aa > 0. In this case we require cr < 0 since cr — 4 < 0 or, because of 

(HZi, 

min (1 - 2^ 

max ■) < -y/Ji max (168) 

which holds if /imax > \ and 

Mmin > M*- (169) 

In this case, for (11781) to hold, a must lie in the range given by (ii). Hence, the 
proof of the theorem is complete. □ 

Corollary 3.5 Consider the GMSSOR method. Let A G and Q G 

j^nxn symmetric positive definite and B G be of full column rank. De¬ 

note the minimum and the maximum eigenvalues of the matrix J = 
by Umin and Umax, respectively. If Umax < \ or if fimax > I and either (i) 
h'min < h* or (a) pimin > and Oi < o < 02 , then the spectral radius of the 
GMSSOR method, p{Q{u!i,ui 2 ,a)), is minimized at 


where 


^opt 


^opt ^2opt 

1 - 


(.y/fi min + y/fi maa:)^ 




(170) 


(171) 


and its corresponding value is 

(nt w y/h-max — y/h-min r7o\ 

P(g(^lopt.^2„pt,a)) = ^ ^ (172) 

Y l-^max “r y j^min 

where pi*, 01 , 02,0 are given by im , [m - 

Proof Recall that GMSSOR is obtained by setting ti = T 2 = w in GMPSD. 
Therefore, (11701) . (I17ip and (|172l) are obtained by (I151L (11521) and (I153p . 
respectively. Indeed, as in GMPSD(3), since = r2„pj it follows that (11771) 
holds also and by the analysis of the proof of Gorollarv l3.41 we have that (I170|) , 
(HZID and (I172|) hold under the same conditions as in Gorollarv l3.4l □ 

Corollary 3.6 Consider the GPSD method. Let A G and Q G 

be symmetric positive definite and B G be of full column rank. Denote 

the minimum and the maximum eigenvalues of the matrix J = Q~^B^A~^B 
by Pmin and pmax, respectively. If pmax < \ or if pimax > J and either (i) 
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Mmin < jJ* or (ii) ^rnin > oud Oi < a < 02, then the speetral radius of the 
GPSD method, p{Q{T,uj,a)), is minimized at 


'^opt — 




and ujopt = 


(\/m min + a/m max ')^ 

and its corresponding value is 

p{Q {t^opt ^opt ^ off) = 

where 


1 ± a/ 1 — a(l — a)a 

\/M max - \/m min 
\/M max + Tm min 


iJ- = 


(l-2y^w)2’ 


Oi = 


= , 02 = 


with 


a = 4(1 — M) and M = 


— 4) (T+a/ cr(cr — 4) 

^P'lnin P'max 


(Vm min + \/m 


(173) 


(174) 


(175) 


(176) 


Proof GPSD follows from GMPSD by letting t = ti = T 2 and uj = oJi = UJ 2 
or Topt = and Uopt = ■ By equating the expressions 

of and r2^pj given by (11521) we obtain 


o(l — o)uJ 2 — UJ 2 + 1 — M = 0, 

where M is given by (ESI). This quadratic has real roots when 

a^a — 0(7 + 1 > 0 , 


(177) 


(178) 


where a is given by (Ei- We distinguish two cases. Gase 1: Aa < 0, Case 2: 
Aa>0 where Aa = (j{a — 4). 

Case 1: Aa < 0. In this case we require cr > 0 since cr — 4 < 0 or in view of 

(fTfSl) 

Vm min (1 - 2^ 

max ') > -\/m max ■ (179) 

But, (11791) holds if either pmax < j or if pmax > \ and Pmin < M* hence (i) is 
proved. 

Case 2: Aa > 0. In this case we require cr < 0 since cr — 4 < 0 or, because of 

ra, 

min (1 - 2^ 

max ■) < -\/m max (180) 

which holds if Pmax > j and 


Mr) 


>M*- 


(181) 


In this case, for (11781) to hold, a must lie in the range given by (ii). Therefore, 
it follows that oJopt is given by (I173|) . □ 

Analogous results hold in case where the matrix Q is symmetric negative 
definite. 










































36 


M. A. Louka et al. 


4 Numerical results 


In this section we study the numerical solution 
equation 

+ Vw = /, 

V-u = 5 , 
u = 0, 
w{x)dx = 0, 

where Q = (0,1) x (0,1) C dfi is the boundary of 17, A is the componen¬ 
twise Laplace operator, u is a vector-valued function representing the velocity 
and ui is a scalar function representing the pressure. Furthermore, we assume 
that the functions f,g are constant. By discretizing (jl82l) with the upwind 
scheme, we obtain the system of linear equations ((H), in which [7] 


of the following linear Stokes 


17 

17 

dn 


(182) 


A = 


iT + T( 
0 


0 

)T-i-r( 


X2p^ 


with 


B = 


[f0iJ 


e 




T = ■ tridiag{—l, 2,-1) S ^ 1,0) S 


h = being the discretization mesh size and 0 the Kronecker product sym¬ 
bol. For this example, we let /i = 1, m = 2p^ and n = . Hence, the total 

number of variables is m -I- n = 3p^. 

We choose the matrix Q to be an approximation to The reason be¬ 

ing that if Q ~ B^A~^B then J = Q~^B^A~^B ~ I. In this case the ratio of 
the maximum to the minimum eigenvalue of the matrix J becomes minimum 
and its value is approximately 1. As a consequence, the spectral radius of the 
iteration matrix of the GMESOR and GMPSD methods attains its minimum 
value. We choose Q, according to the following two cases: 


1. Q = B'^A ^B, A = tridiag(A) 

2. Q = B'^A-^B, A = diag(A), 


where A is the tridiagonal or the diagonal part of A. The choice of the ma¬ 
trix A instead of A is due to the difficulty in computing the inverse matrix of 
A. In this example the eigenvalues of Q are real and positive. In actual com¬ 
putations, we choose the right-hand-side vector (&^, G such that 

the exact solution of the augmented linear system O is Ax*) ,{y*) V = 
(1,1,...,!)^ G and perform all runs in MATLAB (version i?_20126) 

with a machine precision 10“^®. The machine used was an Intel i5 personal 
computer with 6G memory. In our computations, all runs are started from the 
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initial vector 
satisfy 




= 0, and terminated if the current iterations 


V\\b-Ax('^') - By(>=)\\l + \ \q - ^ 

y/\\b- AxW - Byio) \\l + \\q- | |i “ 

where RES is the norm of absolute residual vectors, or if the numbers of the 
prescribed iterations kmax = 1200 are exceeded. We also use the same example 
to compare our methods with the PHSS [7] and Krylov subspace methods m. 

02 ], @9). 

In Table [5] we computed the optimal parameters and and the 

optimal spectral radius popt of the GMESOR method, for various problem 
sizes (m,n) using (ED, (ED and (ED- Furthermore, we computed the opti¬ 
mum parameters T 2 ^^^{exp), uj 2 „pt{exp) and the spectral radii p(T 2 opt{^xp)) 
and p{uj 2 ppt{exp)), experimentally by trial and error. The parameter ti was 
kept fixed and was given its optimum value. Our results show that popt — 
P(^2„pfyea;p)) ~ p{ oj 2 pp,{ exp)) and uj 2 ppt = ~ ~ thus 

verifying Theorem l2.10l The numerical results in Table[7]verify that the param¬ 
eter a may be chosen arbitrary, while the minimum value of p{'H{ti , ra, wa, a)) 
remains approximately the same. , ra„pt, wa„pt, a)) was computed us¬ 

ing Matlab. The slightly different values are due to rounding errors. Finally, in 
Table [5] we list numerical results with respect to the number of total iteration 
steps (denoted by “ITFR”), the elapsed CPU time in seconds (denoted by 
“CPU”) and RES for the CSOR, CMFSOR and Simplified CMPSD iterative 
methods. We remark that our numerical results verify the validity of theorem 
12.101 and corollary 13.31 since CSOR, CMFSOR and Simplified CMPSD meth¬ 
ods require the same number of iterations for convergence. Indeed, this was 
expected since all these methods have the same spectral radius for the opti¬ 
mum values of their parameters. Note that all the aforementioned methods 
require approximately the same computing time. Furthermore, for comparison 
purposes we also considered the PHSS(a*), CMRFS, CMRFS(#), PCMRFS 
and PCMRFS(#) methods. The integer fy in CMRFS(#) and PCMRFS(#) 
methods denotes the number of restarting steps, while the integer a* denotes 
the theoretical optimal parameter of the PHSS method. We also list numerical 
results with respect to the number of total iteration steps and the elapsed CPU 
time in seconds for these methods. The preconditioned matrix Q in PHSS(a*) 
is given by the aforementioned cases 1 and 2. The preconditioner, say K, for 
the PCMRFS and PCMRFS(#) methods is given by [T7], [15], [3^, El 


K = 


A 0 

0 I 


We remark that the CSOR, CMFSOR and Simplihed CMPSD methods al¬ 
ways outperform the other testing methods, except of the PHSS (a*) method, 
considerably with respect to iteration steps as p increases. However, the overall 
computing time of the CSOR, GMESOR and Simplified CMPSD methods is 
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much smaller than that of all the other testing methods. With * we denote 
that the method converges but after too many hours. With regard to the ma¬ 
trix Q, Case 1 is the best choice for all methods tested as it requires the least 
iteration steps and CPU times. 


Table 6 Experimental results showing that u)2opt — GMESOR. 


m 

128 

512 

1152 

n 

64 

256 

576 

m+n 

192 

768 

1728 

Case 1 


0.663309 

0.442911 

0.330674 


0.499375 

0.285422 

0.198468 

‘^2or,t 

0.499375 

0.285422 

0.198468 

popt 

0.580251 

0.746384 

0.811229 

r2o„t [exp) 

0.5 

0.286 

0.199 


0.582936 

0.750508 

0.823517 

[exp) 

0.499 

0.285 

0.198 

p[‘^2ar,t [exp)) 

0.581866 

0.749401 

0.822877 

Case 2 


0.757767 

0.631420 

0.558518 


1.950825 

2.529944 

2.974309 

‘^2or,t 

1.950825 

2.529944 

2.974309 

Popt 

0.492171 

0.607108 

0.664441 

r2o„t [exp) 

1.951 

2.530 

2.975 

MSSKBMStM 

0.492374 

0.607155 

0.664925 

1 ‘^ 2 o„t [exp) 

1.950 

2.529 

2.974 


0.493127 

0.607901 

0.664657 


Table 7 Computation of , T2opt, ^^2opt ’ ^)) various values of the parameter a 

(Case 1, p = 40). 


a 



P['H[Tl^„t , T2„^t . ‘^2o^t > “)) 

0 

2.18851E-001 

1.229935E-001 

0.883807 

10 

2.18851E-001 

5.515564E-002 

0.883808 


2.18851E-001 

9.248083E-003 

0.883808 

■OH 

2.18851E-001 

9.919351E-004 

0.883809 

1 lO'* 

2.18851E-001 

9.991876E-005 

0.883810 


2.18851E-001 

9.999187E-006 

0.883807 

■EH 

2.18851E-001 

9.999919E-007 

0.883807 

1 10' 

2.18851E-001 

9.999992E-008 

0.883807 

■EH 

2.18851E-001 

9.999999E-009 

0.883815 

■EH 

2.18851E-001 

l.OOOOOOE-009 

0.886425 

Hllil 

2.18851E-001 

l.OOOOOOE-010 

0.883879 
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Table 8 ITER, CPU and RES for the testing methods 



P 

8 

16 

24 

32 

40 

48 

n 

64 

256 

576 

1024 

1600 

2304 

m 

128 

512 

1152 

2048 

3200 

4608 

m-|-n 

192 

768 

1728 

3072 

4800 

6912 

Case 1 

GSOR 

ITER 

46 

86 

126 

167 

207 

248 

GPU 

0,05 

0,36 

3771 

22;37 

807 

258,28 

RES 

6.79E-10 

9,04E-10 

9,79E-10 

8,97E-10 

9,74E-10 

9,44E-10 

GMESOR 

ITER 

46 

86 

126 

167 

207 

248 

CPU 

0,05 

0,36 

3,71 

22,61 

86,78 

258,93 

RES 

6.79E-10 

9,04E-10 

9.79E-10 

8,97E-10 

9,74E-10 

9,44E-10 

Simplified 

GMPSD 

ITER 

46 

86 

126 

167 

207 

238 

CPU 

0,05 

0,35 

3,71 

22,59 

86,54 

258,28 

RES 

7,03E-10 

9,12E-10 

9.83E-10 

8,99E-10 

9,75E-10 

9,45E-10 

PHSS(a*) 

ITER 

24 

35 

44 

51 

57 

63 

CPU 

as! 

5123 

33;52 

147,95 

472,87 

1247,80 

RES 

6,19E-10 

9,62E-10 

7.63E-10 

7,36E-10 

9,63E-10 

8,82E-10 

GMRES 

ITER 

73 

176 

285 

386 

506 

606 

CPU 

0,33 

9,24 

155,67 

1.240,62 

6.352,04 

22.214,42 

GMRES(IOO) 

ITER 

73 

327 

831 

1794 

3436 

9965 

CPU 

0,26 

16,12 

404,70 

5.417,99 

41.626,66 

356.831,22 

PGMRES 

ITER 

76 

143 

207 

275 

333 

410 

CPU 

0,50 

11,35 

130,19 

956,58 

4.557,29 

15.684,49 

PGMRES(IOO) 

ITER 

76 

178 

321 

509 

1038 

1281 

CPU 

0,36 

11,02 

172,80 

1.615,68 

12.838,90 

46.595,50 

Case 2 

GSOR 

ITER 

65 

124 

182 

241 

300 

359 

CPU 

ipoT 

032 

4,11 

24,55 

903 

278,67 

RES 

8.35E-10 

8,25E-10 

9,32E-10 

9,14E-10 

9,19E-10 

9,35E-10 

GMESOR 

ITER 

65 

124 

182 

241 

300 

359 

CPU 

0,06 

0,40 

4,08 

24,41 

93,14 

278,48 

RES 

8.35E-10 

8,25E-10 

9,32E-10 

9,14E-10 

9,19E-10 

9,35E-10 

Simplified 

GMPSD 

ITER 

65 

123 

182 

241 

300 

359 

CPU 

0,07 

0,39 

4,11 

24,51 

93,20 

278,08 

RES 

8.55E-10 

8,30E-10 

9,35E-10 

9,15E-10 

9,20E-10 

9,35E-10 

PHSS(a*) 

ITER 

29 

43 

53 

62 

69 

76 

CPU 


5,28 

33;63 

148,51 

474,44 

1261,96 

RES 

9.88E-10 

6,53E-10 

7,99E-10 

8,4SE-10 

9,61E-10 

9,72E-10 

GMRES 

ITER 

73 

176 

285 

386 

506 

606 

CPU 

0,33 

9,24 

155,67 

1.240,62 

6.352,04 

22.214,42 

GMRES(IOO) 

ITER 

73 

327 

831 

r793 

3336 

9965 

CPU 

0,26 

16,12 

404,70 

5.417,99 

41.626,66 

356.831,22 

PGMRES 

ITER 

75 

164 

253 

347 

446 

537 

CPU 

0,50 

12,18 

151,43 

1.168,31 

5.761,78 

20.131,99 

PGMRES(IOO) 

ITER 

75 

275 

544 

997 



CPU 

0,36 

15,31 

279,17 

3.072,70 

> 72h 

» 72h 


5 Remarks and Conclusions 

In this paper we studied the impact of two different preconditioning matrices 
on the convergence of iterative methods for the solution of the augmented lin¬ 
ear system CD when the coefficient matrix A is of the form ©• We assumed 
that A € was a symmetric positive definite matrix and B € 

was a matrix of full column rank, where m > n, whereas Q was a symmetric 
positive or negative definite matrix. Under these assumptions we were able 
to find sufficient conditions for the GMESOR and GMPSD iterative methods 
to converge. Further, using a geometric analysis analogous to Varga m we 
determined the optimum values of the parameters of all methods studied such 
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as to attain the maximum rate of convergence. From our analysis it was shown 
that GMESOR and GMPSD are equivalent since they have the same spectral 
radius for the optimum values of their parameters, which is given by dMl). 
This result was verified by our numerical experiments, where the simplified 
GMPSD, the GMESOR and the GSOR methods require approximately the 
same computing time. Moreover, all the aforementioned methods outperform 
the PHSS(a*), GMRES, GMRES(#), PGMRES and PGMRES(#) methods 
considerably with respect to CPU times. It is worth mentioning that, for the 
saddle point problem, the GMPSD method has a similar behavior as the Mod¬ 
ified PSD (MPSD) method for two-cyclic matrices [33]. Indeed, in [33] we 
proved the equivalence of MPSD and MSOR methods for two-cyclic matrices 
in case the eigenvalues of the Jacobi matrix are either all real or all imaginary. 
However, it is believed that this equivalence will not hold for the case where 
the eigenvalues of the J matrix are complex. 
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